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The PT symmetric potential Vo[cos(27rx/a) + iX sin(27rx/a)] has a completely real 
spectrum for A < 1, and begins to develop complex eigenvalues for A > 1. At the 
symmetry-breaking threshold A = 1 some of the eigenvectors become degenerate, 
giving rise to a Jordan-block structure for each degenerate eigenvector. In general 
this is expected to result in a secular growth in the amplitude of the wave. However, 
it has been shown in a recent paper by Longhi, by numerical simulation and by 
the use of perturbation theory, that for a broad initial wave packet this growth is 
suppressed, and instead a saturation leading to a constant maximum amplitude is 
observed. We revisit this problem by explicitly constructing the Bloch wave-functions 
and the associated Jordan functions and using the method of stationary states to find 
the dependence on the longitudinal distance z for a variety of different initial wave 
packets. This allows us to show in detail how the saturation of the linear growth 
arises from the close connection between the contributions of the Jordan functions 
and those of the neighbouring Bloch waves. 

PACS numbers: 42.25.Bs, 02.30.Gp, 11.30.Er, 42.82.Et 

I. INTRODUCTION 

The study of quantum mechanical Hamiltonians that are PT-symmetric but not 
Hermitian[l|-j6[ has recently found an unexpected application in classical optics0j-[l5[, due 
to the fact that in the paraxial approximation the equation of propagation of an electromag- 
netic wave in a medium is formally identical to the Schrodinger equation, but with different 
interpretations for the symbols appearing therein. The equation of propagation takes the 
form 



where tp(x,z) represents the envelope function of the amplitude of the electric field, z is a 
scaled propagation distance, and V(x) is the optical potential, proportional to the variation 
in the refractive index of the material through which the wave is passing. A complex V 
corresponds to a complex refractive index, whose imaginary part represents either loss or 
gain. In principle the loss and gain regions can be carefully configured so that V is PT 
symmetric, that is V*(x) = V(—x). 
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Propagation through such a medium exhibits many new and interesting properties, such 
as nonreciprocal diffraction [l6[ and birefringence (l3| . One of the main features of complex 
optical lattices is the non-conservation of the total power. In the PT-symmetric case this 
can lead to effects such as power oscillations [l3j. It has been argued that one can distin- 
guish three universal dynamics [TtJ related to broken or unbroken symmetry. While this 
is in general true, the behaviour can be modified considerably for special initial conditions, 
as we will discuss in the present paper. Many familiar effects such as Bloch oscillations 
and dynamical localisation get drastically modified in the presence of imaginary potentials 



and PT-symmetry |18|, [19|. The new features of complex optical lattices provide excit- 
ing opportunities for engineering applications. As an example, the possibility of realizing 
unidirectional light propagation has been envisaged |20| . In the case of high intensities the 
propagation equation ([I]) gets modified due to the Kerr-nonlinearity, leading to an additional 
term proportional to |V ; | 2 ' ? / ; - ft nas been shown in 15] that the influence of the nonlinearity 



on the non-reciprocal effects can be advantageous for applications such as unidirectional 
couplers. It is interesting to note that the nonlinear propagation equation also has a coun- 
terpart in quantum dynamics, as the mean-field description of Bose-Einstein condensates, 



where there has also been interest in PT symmetric models |2l| . However, for the purposes 
of this paper, we shall limit ourselves to the linear case. 

A model system exemplifying some of the novel features of beam propagation in PT- 
symmetric optical lattices uses the sinusoidal potential 

V = Vo [cos(27re/ a) + iX sin(27rx/ a)] . 

This model has been studied numerically and theoretically, e.g. in Refs. 0, [l2[ 13]. The 
propagation in z of the amplitude ip(x, z) is governed by the analogue Schrodinger equation 
(TTJ), which for an eigenstate of H, with eigenvalue (3 and z-dependence ip oc e~ ll3z reduces to 
the eigenvalue equation 

- tp" - V Q [cos(27rx/a) + iX sin(2vrx/a)] ip = /3if; . (2) 

These eigenvalues are real for A < 1, which corresponds to unbroken PT symmetry, where 
the eigenfunctions respect the (anti-linear) symmetry of the Hamiltonian. Above A = 1 
pairs of complex conjugate eigenvalues begin to appear, and indeed above A ~ 1.77687 
all the eigenvalues are complex [221] . Clearly one would expect oscillatory behaviour of the 
amplitude below the threshold at A = 1 and exponential behaviour above the threshold, but 
the precise form of the evolution at A = 1 is less obvious. At first sight one would expect 
linear growth (see, e.g. Ref. because of the appearance of Jordan blocks associated 
with the degenerate eigenvalues that merge at that value of A , but, as Longhi[12j has 
emphasized, this behaviour can be significantly modified depending on the nature of the 
initial wave packet. 

It is this problem that we wish to discuss in the present paper. In Section 2 we explicitly 
construct the Bloch wave-functions and the associated Jordan functions corresponding to 
the degenerate eigenvalues and then use the analogue of the method of stationary states to 
construct the z-dependence. We find that the explicit linear dependence arising from the 
Jordan associated functions is indeed cancelled by the combined contributions from the non- 
degenerate wave-functions (which individually give an oscillatory behaviour). In Section 3 
we analyze this cancellation in detail, showing how the coefficients of the two contributions 
are closely related, and obtaining an approximate analytic expression for the derivative of 
their sum. Our conclusions are given in Section 4. 
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II. BLOCH AND ASSOCIATED JORDAN WAVE-FUNCTIONS 



At the threshold A = 1, the potential V in Eq. (J2j) becomes the complex exponential 
V = Vo exp(2i7rx/a), for which the Schrodinger equation reads 

- V" - V exp(2z7rx/a)^ = /3ip. (3) 

This is a form of the Bessel equation, as can be seen by the substitution y = y exp(i7ix/a), 
where y = (a/n)y/V , giving 

where g 2 = (3(a/ir) 2 . Thus the spectrum is that of a free massive particle, shown in the 
reduced zone scheme in Fig. 1, and for q = ka/ir not an integer the solutions ipk( x ) — Iq{y) 
and ip-k(x) = I- q (y) are linearly independent, and have exactly the correct periodicity, 
ipk{x + a) = e tka ipk{x) , to be the Bloch wave-functions. It is important to note, however, 
that because the original potential is PT-symmetric rather than Hermitian, these functions 
are not orthogonal in the usual sense, but rather with respect to the PT inner product (see 
Eq. ©). 




FIG. 1: Band structure for A = 1 in the reduced zone scheme. The Bloch momentum k is plotted 
in units of 7r/a and the eigenvalue /3 in units of (a/n) 2 . 



A. Jordan Associated Functions 



However, for q = n, a non-zero integer, I n (y) and I- n {y) are no longer independent, 
but are in fact equal, signalling the degeneracy of the eigenvectors at those points, and the 
formation of spectral singularities and Jordan blocks. In that case the Bloch eigenfunctions 
do not form a complete set, and we must search for other functions, still with the same 
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periodicity, to supplement them. These are the Jordan associated functions (see Appendix 
and Refs. [13, 25[ ), which we denote by tpk(x) = Xn{y), defined not by the eigenvalue 
equation itself, but by 



y 



d 2 



d_ 
dy 



Xn(y) = i n (y), 



(5) 



and the periodicity condition Xn{z m y) = e ln Xn(y), corresponding to the Bloch periodicity 
ifik{x + a) = e lka (fk{x). A particular solution of this equation, which can be expressed 
explicitly in terms of generalized hypergeometric functions, is given by 



x F n \y) 



y dz 



I n (z) [K n (z)I n {y) - K n {y)I n {z)\ 



(6) 



as is easily checked by differentiation and use of the Wronskian identity K n (y)I' n (y) — 
K' n (y)I n (y) = 1. However, the corresponding ip^\x) has a discontinuity in its imaginary 
part at x = ±a and does not have the required periodicity. This problem can be rectified by 
recognizing that we may add to Xn 1 an y multiple of I n (y), or more importantly K n (y). The 
latter displays exactly the same kind of discontinuity 1 as x^ 7 , and by choosing its coefficient 
judiciously the discontinuities can be made to cancel. Specifically, we take 



Xn{y) = Xn\y) - 7; K n{y)- 



(7) 



An alternative derivation of this relation using the definition of <fk{x) in terms of d<fik{x)/dk 
will be given in the next section. With the addition of the last term the resulting (fk{x) is 
not only free of discontinuities, but also obeys the correct Bloch periodicity condition. This 
is shown in Figure 2, where we plot the real and imaginary parts of y9^ 7 (x) and tpk( x ) for 
q = 2. Note the PT symmetry of ipk(x), namely y?fc(— x) = fl(x). 





(a) 

FIG. 2: (color online) The real (green, continuous) and imaginary (red, dashed) parts of (a) the 
particular integral ip? 1 \x) and (b) the corrected Jordan function (fk{%) for a = ir,k = q = 2 and 
V = 2. 



This arises because K n (y) has a branch cut along the negative real axis, and when x passes through ±a 
one should really go to the next Riemann sheet. 
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The complete set of functions, orthogonal with respect to the PT metric, now consists of 
the Bloch eigenf unctions ipk( x ), supplemented by the Jordan associated functions <p k (x) for 
k = rnr/a with n > 0, and a general wave-function f(x) may be expanded as 

f(x) = c i>o(x) + ^ c ktpk(x) + ^2(a n I n (y) + n Xn(v))- (8) 

kj^mr/a n> 

Here we have discretized the problem by putting the system in a box of length 2Na, in 
which case k — )■ k r = rir/(Na). The coefficients c k can be obtained[9j by using the PT 
orthogonality 



dxi/}- k (x)ip k ,(x) = 5 kk , / dxip_ k (x)ip k (x), (9) 



as 



Ck = j dxj)- k {x)f{x) 
k J dxip_ k (x)ip k {x) 

Here the sign of the denominator alternates from band to band, reflecting the indefinite 
nature of the PT metric. However, this relation breaks down precisely at the Brillouin 
Zone boundaries, where the single Bloch eigenfunctions are self-orthogonal, J dx I^iy) = 0, 
another indication that we need the supplementary Jordan functions. 

Thus it is P n , rather than a n , which is determined by integrating f(x) with respect to 
l n {y): 

p = Jdxl n {y)f{x) 
f dxI n (y)xn(yY 

while a n is subsequently determined by integrating f(x) with respect to Xniv)- 
As an immediate check of the correctness of the functions Xm th e identity 

1 = I (y) - 21 2 {y) + 2I 4 (y) (12) 

implies that 

J dx X 2n{y) = 2(-l) n J dxI 2n (y) X 2n(y), (13) 
for n > 0, a relation we have verified numerically. 



B. Method of Stationary States 

In standard quantum mechanics one method of solving the time-dependent Schrodinger 
equation 

tj t ^(x,t) = H^(x,t) (14) 

is to expand the initial wave-function ip(x, 0) in terms of the (complete) set of orthonormal 
eigenfunctions ip r as 

ip(x, 0) = a m tp m (x), (15) 
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with the coefficients a m given by the overlap 

dx^^ n (x)ij;(x,0). (16) 

The wave-function at time t is then given by 

V(x,t) = ^a m e-^Vm(x). (17) 

m 

Here we have essentially the same problem, with z taking the role of t, but with the crucial 
difference that we must include the Jordan associated functions in the sum, which now takes 
the form of Eq. dHJ), with the coefficients determined as in Eqs. (ITO]) and (TTT]) . As is well 
known, and has been emphasized again recently in Refs. [12| and|23||. the time dependence 
then takes a different form. The Jordan function <p r satisfies the equation (H — E r )<p r = ip r , 
and hence 

= e- iErt {ip r -it^ r ). (18) 

Thus, in addition to the usual phase factors making up the sum for the time-dependent wave- 
function, one has an explicit factor of t multiplying the degenerate eigenfunctions associated 
with a Jordan block, giving the complete time dependence of the initial wave-function of 
Eq. © as 

-ik 2 t 



f{x,t) = c ip (x)+ ^ Ck^k{x)e 

kj^nn/a 

+ £ (K - it(3 n )I n {y) + (3 nXn (y))e-^ 2/a2)t . (19) 



kj^nn/a 

„2„2 i„2\ 



n> 



However, the explicit factor of t only appears when the coefficient of the Jordan function is 
nonzero. The expansion of Eq. (j!2p is a case in point. 

In the corresponding optical problem it therefore appears that the z-dependence should 
be expected to be oscillatory for A < 1, exponential for A > 1 and linear precisely at the 
symmetry-breaking threshold A = 1 for initial states that excite a Jordan function. However, 
Longhi Jjl] has shown numerically and in perturbation theory that this linear dependence 



may not be realized, depending on the nature of the initial wave-function. Since we now have 
explicit expressions for the Bloch and Jordan functions we are in a position to investigate 
the origin of this phenomenon in the context of the method of stationary states. 
We will take as our input a Gaussian profile of the form 

iP(x, 0) = f(x) = e - {x/w)2+ikoX , (20) 

with offset fco and width w. Because of the periodicity of the Bloch eigenfunctions and 
Jordan associated functions the range of the integral in Eqs. ffTOl) and ffTTl) can be reduced 
to < x < a, provided that f(x) is replaced by 

N-l 

F q (x)= J2 e- inmq f{x + ma), (21) 

m=-N 



7 



in which we recall that q = ka/n. That is, 



/° dxip- k (x)Fg(x) 
° k 2N ft daaJ>- k (x)M*V 
and similarly for a n , (3 n . In particular 

£dxI n (y)F n (x) 

Pn 2N j a dxiMxM' [ } 

For N large, the sum in Eq. (121 j) can be extended to infinity without significant error. 
However, we have to be careful about the periodicity in q of the discretized version 2 . That 
is, F q as given in Eq. (j2ip is really a function of g mod 2. Thus we first replace q by g, 
where q is the nearest to go of the set of equivalent momenta, satisfying |g — go| < 1, to 
obtain the result 

F q (x) = e -i- 2 (^o) W+-^ 3 L X _ + l i;r 2 (g - _ ^ j e -^>A ; (24) 



where $3(2:, v) is the Jacobi elliptic theta function, which has the expansion 26 



^ 3 (z,v) = 1 + 2^2 v s2 cos(2sz). (25) 



2 2/2 

For a broad wave-packet, with w ^> a, the argument t> = e _7r ™ > a is very small, so that 
$3(2, ~ 1. Moreover the prefactor in Eq. (j24p is also small except for the case q = go = 
fcoa/TT- 

We are now in a position to identify under what conditions the Jordan associated functions 
are excited, that is, when /3 n , as given by Eq. f[2Uj) . is non- vanishing. In order to obtain an 
appreciable value of F n the scaled offset momentum g must be an integer m, with n = m 
mod 2. Now, however, there is the integral in the numerator to be considered. Given that 
the phase of F n (x) is approximately e mqx ^ a and that I n (y) has the expansion 

1 00 (±v 2 ) s 
W = M)°E s , r( l 4 !I +1) . P6) 



=0 



where we recall that y = yoe lnx ^ a , it is easily seen that the integral vanishes unless m is 
negative and n < \m\. In that case the Jordan functions <p\ m \, f\m\-2, ■■■ are excited. 
However, if go is not a negative integer, no Jordan blocks are excited, and thus no linear 
growth is to be expected. For g = no Jordan function is excited because the ground- 
state level n = is non-degenerate. Thus no linear growth is expected in this case. Nor 



is it expected for the case g = +1. In the set-up chosen by Longhi in Ref. [121 ]. on the 
other hand, the offset go was taken to be -1, in which case ip\ is excited. Note the left-right 
asymmetry here: the Hamiltonian is not parity invariant, but only PT invariant. 

Figures 3 and 4 show the different propagation behaviour for a wide beam in the two 
cases go = — 1 and go = respectively. The parameters are: a = 11, w = 6ir and Vq = 2. 



2 This is similar to the difference between the delta function and its discretized version, the sine function. 
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In the first case, where one Jordan mode tpi, is excited, the beam spreads out but does not 
split, as shown in Fig. 3(a). In Fig. 3(b) we show the different contributions to the maximum 
amplitude. The lower curve shows the contribution of the Jordan block sector only, which 
indeed rises linearly, as expected. However, the intermediate curve, the contribution of all 
other modes, mysteriously begins to decrease after an initial rise, and the upper curve, which 



takes into account all the contributions, exhibits the saturation first noted by Longhi|l2 

(a) 





FIG. 3: (color online) (a) \i(i{x, z)\ as a function of x, z. (b) Maximum value (red solid line) of 
|?/>(x,z)| as a function of z. Blue dashed-dotted line: Jordan block contributions (second line of 
Eq. (|19p only). Green dashed line: other contributions only. The parameters are: a = n, Vq = 2, 
w = 6-7T and qo = —1. 




In the second case, where no Jordan mode is excited, the behaviour is quite different. The 
beam does not significantly spread, but instead splits into two, exhibiting the phenomenon 
of birefringence [9|, and the maximum amplitude, after some initial fluctuations, decreases 
slowly with z. The behaviour for the case when go — 1 is similar. 



III. THE MECHANISM OF SATURATION 



In this section we investigate in detail how the contributions from the other (non- Jordan) 
eigenfunctions conspire to cancel out the explicit linear growth in z coming from the Jordan 
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sector. A priori such a cancellation seems highly unlikely, but since it does occur it must 
be because the two sets of contributions are in fact closely connected. Mathematically 
the cancellation cannot be complete, but can only happen only over a limited range in z. 
However, that is sufficient for the physical situation, because our lattice is finite in the x 
direction, and at very large values of z the beam will have encountered the edges of the 
lattice. 

We now show that the contributions of the Jordan block and of the nearby eigenfunctions 
are indeed closely related. We will consider specifically the first case of the previous section, 
where the Jordan mode cpx is excited. 

Recall that ip(x, 0) = f(x) is expanded as 

lf)(x, 0) = C l/J (x) + ^ °k^k(x) + ^2 ( a n!n(y) + &nXn{y)), (27) 
k^nn/a n> 

as in Eq. ([S]). The expression for c k is given in Eq. (fIU|) . in which the integration over x 
can be reduced to the interval [0, a] by exploiting the periodicity of the Bloch functions, to 
obtain 

^dx'ip_ k (x)F q (x) 

Ck — nl . T r a j ; , \ , / \ i \ z °) 



2N J Q a dxip- k (x)ip h (x) 



where 



N-l 



F q (x)= e- inmq f(x + ma). (29) 



m=-N 



Similarly (3 n is given by 



/ a dxI n (y)F n (x) 

Pn 2N£ dxl n (y)xn(y) ' [ ' 

The denominator of would vanish at q = n. In fact it turns out that it is precisely a 
sine function: 

(■a 

dxip-k{x)ipk{x) = a sinc(ka) = a sinc(g7r). (31) 

The denominator of (3 n is proportional to the derivative of this previous denominator with 
respect to k. Thus 

d f a a d f a 

— / dxifj_ k (x)ip k (x) = -— / dxl- q {y)lq{y) 
dk ./n vr da ./ n 



[ a dx (lL q (y)I q (y) + l- q (y)W) ■ ( 32 ) 

Jo 



But by differentiating the general relation I- q (y) = I q {y) + (2/7r) sin qn K g (y) and setting 
q = 1 we find that ILi(y) = I[(y) — Ki{y). Thus, using the derivative definition of the 
Jordan associated function as 3 y?f 7 = (1/ (2k))dipk/dk, we see by reference to Eq. ([7]) that 



Ik I (/ - r ' '-/-■'>')' '/, <>•) 



^ = ^J*dxh(y) (*f'(v) - i#i(v)) 



3 In fact this differs from ip PI as given in Eq. © by a multiple of I q (y) 
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= 4- / dxh(y) X i(y) ■ (33) 

JO 

The numerator of j3 n is a smooth continuation of the numerator of c k , which from now on 
we denote as cV Moreover, is highly peaked around q = 1 (and also around q — — 1), as 
is shown in Fig. 5. In fact, near g = 1 it is given by the Gaussian dk oc e _7r e w /( 4a ) ; where 
e = q — 1, while the denominator is proportional to e. 

A 

c k 



0.30 


• 


0.25 




0.20 




0.15 




0.10 




0.05 




£ 1 l--.--, 


• • 



FIG. 5: The numerator c/% of the expansion coefficients c& in Eq. (|28|) for ^ = 40, a = tt, w = 6tt 
and Vo = 2. At the points q = ±1, which are not included in the expansion, the figure gives instead 
the value of the numerator of f3\ 



We are now in a position to examine the ^-development of the Jordan contribution and 
the related contributions from nearby values of q. Recall that ip(x, z) is given by 

ip(x,z) = c ^o(x)+ ^2 c k i/j k (x)e~ lh2z 

kj^nn/a 

+ IK - Wn)Uy) + PnXn{y)]e-^ 2/a2)X . 
n> 

Thus the total contribution from the neighbourhood of q = 1 and q = — 1 is 

ip(x,z)j ps const x h(y)e 

where e r = r/N. In principle the upper limit of the sum is infinity, but in practice it can be 
taken less than N. Note that if e were continuous the limit of the second term as e — > would 
be twice the first term, which comes from the Jordan function. It turns out that although 
\ip(x,z)j\ does exhibit the expected linear behaviour in z initially, it subsequently|-?/>(a;, z)j\ 
has an extremely wide and flat plateau before it eventually rises again as z approaches Ntt. 
This is illustrated in Fig. 6(a) for N = 40. A hint of such a plateau- like behaviour can be 
seen in the simple function z + sin(2z)/2, which is plotted in Fig. 6(b). However, in this 
case the plateau is much less pronounced. 



ir 2 e r 



sin 



2tt 2 



-e r z 



-el-K 2 (w 2 /A+iz)/a 2 



(34) 
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(a) (b) 

FIG. 6: (a) Modulus of the function in square brackets in Eq. (|34p as a function of z for N = 40, 
a = 7r and w = 6ir. (b) The function z + sin(2z)/2. 



We can understand the extreme flatness of the plateau in Fig. 6(a) in terms of Jacobi $3 
functions. For simplicity let us set a = 7r, so that we are analyzing the function 



+ J2 (-) sin(2e r z)e- e '/^ 2 /4+te) (35) 

r=l \ €r ' 

+ Yl ( ~~) Sin ( 2e r 
r=l ^ r ' 



00 



for the values of z we are considering. While this is not itself a $3 function, its derivative 
with respect to z can be so expressed: 

00 

f(z) « l + 2^cos(2e r z)e- e ' u ' 2/4 



r=l 

00 



1 + 2^COs(2r,^/iV)e- 



■r 2 u; 2 /( 4A^2 ) 



r=l 



= ^^,e- w2/(4JV2) J • (36) 

The behaviour of this $ function is not immediately apparent, since the second argument 
is of 0(1) for w <C 2N. However, it can be made clear by using the alternative notation 
$(z,q) = "&s(z\t), where q = e mr , and applying Jacobi's imaginary transformation 27|. 
whereby 

■& z {z\t) = (-ir)-h- iT ' z2/{nT,) M^'\r'), (37) 
where r' = — 1/r. This converts j'{z) to 

j\z) = 2 v ^-e-^/-^ 3 f e-^ N2 'A . (38) 



In this form the second argument of $3 is small, so that for moderate values of z we can 
approximate $3 by 1, in which case the behaviour is dominated by the preceding Gaussian, 
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which rapidly falls from 1 to a very small value, corresponding to the plateau in j(z). The 
Gaussian is eventually overwhelmed by the hyperbolic cosines occurring in the expansion of 
$3, as must be the case, since $3 is periodic in z. The plot of this function is given in Fig. 7. 

6 

4 
2 



20 40 60 80 100 120 
FIG. 7: j'(z) from Eq. ([38]) versus z. The parameters are the same as in Fig. 6. 



IV. CONCLUSIONS 

On general grounds one expects linear growth of the amplitude ip(x, z) at the symmetry- 
breaking threshold A = 1, due to the degeneracy of a subset of the eigenf unctions at this 
point and the consequent development of Jordan blocks. However, it has been observed 
numerically (l2| that this growth becomes saturated at large z, at least for wide input beams. 
We have been able to explain this saturation phenomenon by analyzing in detail the separate 
contributions from the Jordan blocks and the contributions from the nearby Bloch functions 
in which they are embedded. 

For the particular potential considered here, the Bloch eigenfunctions are associated 
Bessel functions of the first kind, I q (y), and we were able to explicitly construct the associated 
Jordan functions at the exceptional points q = n. Hence we were in a position to use the 
analogue of the method of stationary states to generate the z dependence. In this method 
we were able to isolate the separate contributions of the Jordan blocks from the other non- 
degenerate states with normal z dependence generated by multiplying each eigenfunction 
by its appropriate z-dependent phase. We found that, in cases where the associated Jordan 
functions are excited (Fig. 3) the explicit linear increase of the Jordan-block contributions is 
precisely compensated by a linear decrease of the contribution of the non-degenerate states, 
which of course have no explicit linear ^-dependence. In cases where the Jordan associated 
functions are not excited (Fig. 4), there is no initial linear increase, but rather a rapid 
oscillatory behaviour followed by a very slow decrease. The topology of the beams in these 
two cases is markedly different. In the first case, although the maximum amplitude becomes 
constant, the beam spreads laterally in a linear fashion, and hence the total power grows 
linearly, with the beam taking energy from the lattice. In the second case the total power 
turns out to be constant for large z, with the slow decrease in the maximum amplitude being 
matched by the slow broadening of the individual beams. 
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We then examined in detail the mechanism of saturation, which is only possible because of 
the close relation between the contribution of the Jordan block and the contributions of the 
Bloch functions in which it is embedded. With the aid of a certain amount of approximation 
we were able to write the respective contributions in the relatively simple form of Eq. ( I34p 
and to express its z-derivative as a $3 function. This $3 function encodes the initial increase 
and the subsequent extremely wide and flat plateau. As a mathematical function it also 
encodes further steps and plateaux for larger z, but these are not physically relevant because 
they correspond to values of x outside the range of the finite lattice. 



Acknowledgments 

We are grateful to Prof. R. J. Rivers for extremely useful comments. EMG is supported 
by an Imperial College Junior Research Fellowship. 



APPENDIX 

Since the references for Jordan associated functions are not readily accessible, we give 
here a brief outline of their main properties. They are the analogue of the Jordan associated 
vectors that occur in simple matrix eigenvalue problems of the form 

(H - X)u = (39) 

where H is non-Hermitian. For Hermitian problems, while two eigenvalues Ai and A2 may 
become degenerate at a particular critical value of a parameter in the Hamiltonian, there 
remain two distinct eigenvectors. However, in the non-Hermitian case it is possible that the 
eigenvectors also become degenerate: U\ — U2- The reduction in the number of eigenvectors 
at such a point means that the set of eigenvectors no longer forms a complete basis. However, 
the basis can be completed by inclusion of the Jordan associated vector v\ (generalized 
eigenfunction), defined by the generalized eigenvalue equation 

{H - X 1 )v 1 = Ul . (40) 

Of course this associated vector is not uniquely defined: to it may be added any multiple of 

Ui. 

The archetypical example is given by the non-Hermitian Jordan matrix 

m =(oa)' (41) 

which has the single eigenvalue A and only one eigenvector u = (1,0). The independent 
vector v = (0, 1) needed to complete the basis is indeed a solution (undetermined up to a 
multiple of u) of the generalized eigenvalue equation 

(M — X)v = u . (42) 

The eigenvectors of a non-Hermitian matrix are not orthogonal in the usual sense. Instead 
one needs to introduce the left eigenvectors Ul, satisfying Ul(H — A) = 0, which are different 
from the usual (right) column vectors satisfying (H — \)ur = 0. The orthogonality is then 



14 



between left and right eigenvectors: (ul)i(ur)^ = for distinct eigenvalues Ai and A 2 . 
When the eigenvectors become degenerate, they are self-orthogonal. This is exemplified by 
the Jordan matrix M in Eq. (|4ip . whose left and right eigenvectors are ul = (0,1) and 
ur = u= (1,0), respectively, which are indeed orthogonal. 

It is easy to show in general that the associated vectors v n are orthogonal in the same 
sense to eigenvectors u m and associated vectors v m belonging to different eigenvalues. In 
the event that H has a PT symmetry, for some definition of the reflection operator P, all 
of these overlaps can instead be expressed in terms of right eigenvectors alone with the aid 
of the PT metric. 

The continuum analogue of this situation is the eigenvalue problem 

(H - E n )^j n = (43) 

where H is a non-Hermitian differential operator with eigenvalue E n and (right) eigenfunc- 
tion ip n . Again, it is possible that, at a particular critical value of a parameter in H two 
eigenvalues E m and E n , and also the corresponding eigenfunctions ip m and ip n , become de- 
generate. At that point the single eigenfunction becomes self-orthogonal with respect to the 
metric defined by H, which, for a PT-symmetric problem, is the PT metric. 

The reduction in the number of eigenfunctions at such a point again means that the 
set of eigenfunctions no longer forms a complete basis. However, in complete analogy with 
Eq. (f42|) . the basis can be completed by inclusion of the Jordan associated functions ip n 
(generalized eigenfunctions), defined by the generalized eigenvalue equation 

(H - E n )cp n = ij; n . (44) 

This is precisely how the Jordan functions were introduced in Eq. They are defined 
only up to multiples of solutions of the homogeneous equation and are guaranteed to be 
orthogonal, using the PT metric, to each other and to eigenfunctions belonging to other 
eigenvalues. We found it necessary to exploit this freedom in Eq. ([7]) in order to ensure that 
the ip n satisfy the appropriate boundary conditions. 

An alternative definition of the Jordan functions is as derivatives of the eigenfunctions 
with respect to the energy [2^]. Thus, differentiating the eigenvalue equation (H — E)ip = 
with respect to E we get precisely 

{H-E)^-1> = 0, (45) 

so that we may identify (p with dip/dE, again modulo solutions of the homogeneous equation. 

In our present problem this latter definition leads to extremely simple formulae for the 
functions \ n . For example, using Eq. (9.6.44) of Ref. [26[ it yields Xi{u) — — /o(2/) / (2?/) , 
which is easily seen to be a solution of Eq. ([5]) with the correct periodicity. 
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